﻿/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#ifndef _GCTL_POINT3S_H
#define _GCTL_POINT3S_H

#include "../core/macro.h"
#include "../core/exceptions.h"
#include "point3c.h"

#include "iostream"
#include "string"
#include "cmath"
#include "iomanip"
#include "regex"

namespace gctl
{
	template <typename T> struct point3c;
	template <typename T> struct point3s;

	typedef point3s<double> point3ds;
	typedef point3s<float>  point3fs;

#ifndef IO_PSN
#define IO_PSN
	// static variable for controlling the IO process
	static int io_psn = 6;
#endif // IO_PSN

	/**
	 * @brief      球坐标系中的一个点（经纬度表示的）
	 */
	template <typename T>
	struct point3s
	{
		T rad, lon, lat;
		// 构造函数与拷贝构造函数
		point3s();
		point3s(T r, T ln, T lt);
		point3s(const point3s<T> &b);
		// 析构函数
		virtual ~point3s(){}
		// 结构体方法
		bool valid() const;
		void set(T r, T ln, T lt);
		void set(const point3s<T> &b);
		double module() const;
		point3c<T> s2c() const; // 转换为球坐标点
		void str(std::string str);
		/**
		 * @brief      Sets the i/o precision of the type.
		 *
		 * @param[in]  psn   The desired precision
		 */
		void set_io_precision(int psn);
		/**
		 * @brief      输出位置
		 *
		 * @param      os    输出流
		 */
		void out_loc(std::ostream &os, char deli) const;
		/**
		 * @brief      输入位置
		 *
		 * @param      os    输入流
		 */
		void in_loc(std::istream &os);
		/**
		 * @brief      返回位置
		 *
		 * @return     位置
		 */
		point3s<T> get_loc() const;
		/**
		 * @brief 返回正常球坐标中的Phi值
		 * 
		 * @return double 
		 */
		T get_phi() const;
		/**
		 * @brief 返回正常球坐标中的Theta值
		 * 
		 * @return double 
		 */
		T get_theta() const;
	};

	template <typename T>
	gctl::point3s<T>::point3s()
	{
		rad = lon = lat = NAN;
	}

	template <typename T>
	gctl::point3s<T>::point3s(T r, T ln, T lt)
	{
		set(r, ln, lt);
	}

	template <typename T>
	gctl::point3s<T>::point3s(const point3s<T> &b)
	{
		set(b);
	}

	template <typename T>
	bool gctl::point3s<T>::valid() const
	{
		if (std::isnan(rad) || std::isnan(lon) || std::isnan(lat)) return false;
		if (std::isinf(rad) || std::isinf(lon) || std::isinf(lat)) return false;
		if (rad < 0 || lon > 180 || lon < -180 || lat > 90 || lat < -90) return false;
		return true;
	}

	template <typename T>
	void gctl::point3s<T>::set(T r, T ln, T lt)
	{
		if (std::isnan(r) || std::isnan(ln) || std::isnan(lt) || 
			std::isinf(r) || std::isinf(ln) || std::isinf(lt))
		{
			throw invalid_argument("Invalid value detected. From point3ds::set(...)");
		}

		if (r < 0 || ln > 180 || ln < -180 || lt > 90 || lt < -90)
		{
			throw out_of_range("Invalid values. From point3s::set(...)");
		}

		rad = r; lon = ln; lat = lt;
		return;
	}

	template <typename T>
	void gctl::point3s<T>::set(const point3s<T> &b)
	{
		if (std::isnan(b.rad) || std::isnan(b.lon) || std::isnan(b.lat) || 
			std::isinf(b.rad) || std::isinf(b.lon) || std::isinf(b.lat))
		{
			throw invalid_argument("Invalid value detected. From point3s::set(...)");
		}

		if (b.rad < 0 || b.lon > 180 || b.lon < -180 || b.lat > 90 || b.lat < -90)
		{
			throw out_of_range("Invalid radius. From point3s::set(...)");
		}

		rad = b.rad; lon = b.lon; lat = b.lat;
		return;
	}

	template <typename T>
	double gctl::point3s<T>::module() const
	{
		return rad;
	}

	template <typename T>
	gctl::point3c<T> gctl::point3s<T>::s2c() const
	{
		/*
		point3c<T> outc;
		outc.x = rad*sin((0.5 - lat/180.0)*GCTL_Pi)*cos((2.0 + lon/180.0)*GCTL_Pi);
		outc.y = rad*sin((0.5 - lat/180.0)*GCTL_Pi)*sin((2.0 + lon/180.0)*GCTL_Pi);
		outc.z = rad*cos((0.5 - lat/180.0)*GCTL_Pi);
		return outc;
		*/

		point3c<T> v;
		v.x = rad*cos(GCTL_Pi*lat/180.0)*cos(GCTL_Pi*lon/180.0);
		v.y = rad*cos(GCTL_Pi*lat/180.0)*sin(GCTL_Pi*lon/180.0);
		v.z = rad*sin(GCTL_Pi*lat/180.0);
		return v;
	}

	template <typename T>
	void gctl::point3s<T>::str(std::string str)
	{
		std::smatch ret;
		std::regex pattern("\\((-?\\d*\\.?\\d+?),[ ]*(-?\\d*\\.?\\d+?),[ ]*(-?\\d*\\.?\\d+?)\\)");

		if (regex_search(str, ret, pattern))
		{
			rad = atof(std::string(ret[1]).c_str());
			lon = atof(std::string(ret[2]).c_str());
			lat = atof(std::string(ret[3]).c_str());
			return;
		}

		throw runtime_error("Fail to parse the input string: " + str + ". From point3s::str(...)");
	}

	template <typename T>
	void gctl::point3s<T>::set_io_precision(int psn)
	{
		if (psn < 0)
		{
			throw invalid_argument("Invalid precision. From point3s::set_io_precision(...)");
		}

		io_psn = psn;
		return;
	}

	template <typename T>
	void gctl::point3s<T>::out_loc(std::ostream &os, char deli) const
	{
		os << std::setprecision(io_psn) << rad << deli << lon << deli << lat;
		return;
	}

	template <typename T>
	void gctl::point3s<T>::in_loc(std::istream &os)
	{
		os >> rad >> lon >> lat;
		return;
	}

	template <typename T>
	gctl::point3s<T> gctl::point3s<T>::get_loc() const
	{
		return point3s<T>(rad, lon, lat);
	}

	template <typename T>
	T gctl::point3s<T>::get_phi() const
	{
		if (lon < .0) return GCTL_Pi*(2.0 - lon/180.0);
		return GCTL_Pi*lon/180.0; 
	}

	template <typename T>
	T gctl::point3s<T>::get_theta() const
	{
		return GCTL_Pi*(0.5 - lat/180.0);
	}

	template <typename T>
	bool operator ==(const point3s<T> &a, const point3s<T> &b)
	{
		if(fabs(a.lon-b.lon)<=GCTL_ZERO && fabs(a.lat-b.lat)<=GCTL_ZERO && fabs(a.rad-b.rad)<=GCTL_ZERO)
		{
			return true;
		}
		else return false;
	}

	template <typename T>
	bool operator !=(const point3s<T> &a, const point3s<T> &b)
	{
		if(fabs(a.lon-b.lon)>GCTL_ZERO || fabs(a.lat-b.lat)>GCTL_ZERO || fabs(a.rad-b.rad)>GCTL_ZERO)
		{
			return true;
		}
		else return false;
	}

	template <typename T>
	point3s<T> operator -(const point3s<T> &a, const point3s<T> &b)
	{
		point3s<T> m;
		m.rad = a.rad - b.rad;
		m.lon = a.lon - b.lon;
		m.lat = a.lat - b.lat;

		while (m.lon > 180) m.lon -= 360;
		while (m.lon < -180) m.lon += 360;
		while (m.lat > 90) m.lat = 180 - m.lat;
		while (m.lat < -90) m.lat = -180 - m.lat; 
		return m;
	}

	template <typename T>
	point3s<T> operator +(const point3s<T> &a, const point3s<T> &b)
	{
		gctl::point3s<T> m;
		m.rad = a.rad + b.rad;
		m.lon = a.lon + b.lon;
		m.lat = a.lat + b.lat;

		while (m.lon > 180) m.lon -= 360;
		while (m.lon < -180) m.lon += 360;
		while (m.lat > 90) m.lat = 180 - m.lat;
		while (m.lat < -90) m.lat = -180 - m.lat; 
		return m;
	}

	template <typename T>
	point3s<T> operator *(int sign, const point3s<T> &b)
	{
		gctl::point3s<T> m;
		m.rad = sign*b.rad;
		m.lon = sign*b.lon;
		m.lat = sign*b.lat;

		while (m.lon > 180) m.lon -= 360;
		while (m.lon < -180) m.lon += 360;
		while (m.lat > 90) m.lat = 180 - m.lat;
		while (m.lat < -90) m.lat = -180 - m.lat; 
		return m;
	}

	template <typename T>
	point3s<T> operator *(float sign, const point3s<T> &b)
	{
		gctl::point3s<T> m;
		m.rad = sign*b.rad;
		m.lon = sign*b.lon;
		m.lat = sign*b.lat;

		while (m.lon > 180) m.lon -= 360;
		while (m.lon < -180) m.lon += 360;
		while (m.lat > 90) m.lat = 180 - m.lat;
		while (m.lat < -90) m.lat = -180 - m.lat; 
		return m;
	}

	template <typename T>
	point3s<T> operator *(double sign, const point3s<T> &b)
	{
		gctl::point3s<T> m;
		m.rad = sign*b.rad;
		m.lon = sign*b.lon;
		m.lat = sign*b.lat;

		while (m.lon > 180) m.lon -= 360;
		while (m.lon < -180) m.lon += 360;
		while (m.lat > 90) m.lat = 180 - m.lat;
		while (m.lat < -90) m.lat = -180 - m.lat; 
		return m;
	}

	template <typename T>
	point3s<T> operator *(const point3s<T> &b, int sign)
	{
		gctl::point3s<T> m;
		m.rad = sign*b.rad;
		m.lon = sign*b.lon;
		m.lat = sign*b.lat;

		while (m.lon > 180) m.lon -= 360;
		while (m.lon < -180) m.lon += 360;
		while (m.lat > 90) m.lat = 180 - m.lat;
		while (m.lat < -90) m.lat = -180 - m.lat; 
		return m;
	}

	template <typename T>
	point3s<T> operator *(const point3s<T> &b, float sign)
	{
		gctl::point3s<T> m;
		m.rad = sign*b.rad;
		m.lon = sign*b.lon;
		m.lat = sign*b.lat;

		while (m.lon > 180) m.lon -= 360;
		while (m.lon < -180) m.lon += 360;
		while (m.lat > 90) m.lat = 180 - m.lat;
		while (m.lat < -90) m.lat = -180 - m.lat; 
		return m;
	}

	template <typename T>
	point3s<T> operator *(const point3s<T> &b, double sign)
	{
		point3s<T> m;
		m.rad = sign*b.rad;
		m.lon = sign*b.lon;
		m.lat = sign*b.lat;

		while (m.lon > 180) m.lon -= 360;
		while (m.lon < -180) m.lon += 360;
		while (m.lat > 90) m.lat = 180 - m.lat;
		while (m.lat < -90) m.lat = -180 - m.lat; 
		return m;
	}

	template <typename T>
	std::ostream &operator <<(std::ostream & os, const point3s<T> &a)
	{
		os << std::setprecision(io_psn) << a.rad << " " << a.lon << " " << a.lat;
		return os;
	}

	template <typename T>
	std::istream &operator >>(std::istream & os, point3s<T> &a)
	{
		os >> a.rad >> a.lon >> a.lat;
		return os;
	}

	/**
	 * @brief      函数判断两个point3ds类型是否等于
	 *
	 * @param[in]  a     三维空间内的一个向量或实点a。
	 * @param[in]  b     三维空间内的一个向量或实点b。
	 * @param[in]  cut_off  截断误差，默认值为 gctl_macro.h 中预设的 GCTL_ZERO 值，可在调用时改变。
	 *
	 * @return     是否相等。
	 */
	template <typename T>
	bool isequal(const point3s<T> &a, const point3s<T> &b, double cut_off = GCTL_ZERO)
	{
		if(fabs(a.rad-b.rad) <= cut_off && fabs(a.lon-b.lon) <= cut_off && fabs(a.lat-b.lat) <= cut_off)
		{
			return true;
		}
		else return false;
	}

	/**
	 * @brief      函数判断两个point3ds类型是否不等于
	 *
	 * @param[in]  a     三维空间内的一个向量或实点a。
	 * @param[in]  b     三维空间内的一个向量或实点b。
	 * @param[in]  cut_off  截断误差，默认值为 gctl_macro.h 中预设的 GCTL_ZERO 值，可在调用时改变。
	 *
	 * @return     是否不相等。
	 */
	template <typename T>
	bool notequal(const point3s<T> &a, const point3s<T> &b, double cut_off = GCTL_ZERO)
	{
		if(fabs(a.rad-b.rad)>cut_off || fabs(a.lon-b.lon)>cut_off || fabs(a.lat-b.lat)>cut_off)
		{
			return true;
		}
		else return false;
	}

	/**
	 * @brief      生成直球坐标系的点数组
	 *
	 * @param      out_ps  返回的点数组
	 * @param[in]  lonmin  lon最小值
	 * @param[in]  lonmax  lon最大值
	 * @param[in]  latmin  lat最小值
	 * @param[in]  latmax  lat最大值
	 * @param[in]  dlon    lon间隔
	 * @param[in]  dlat    lat间隔
	 * @param[in]  rad     半径值
	 */
	template <typename T>
	void get_grid_point3s(array<point3s<T>> &out_ps, T lonmin, T lonmax, T latmin, 
		T latmax, T dlon, T dlat, T rad)
	{
		if (lonmin >= lonmax || latmin >= latmax || lonmin+dlon>lonmax || latmin+dlat>latmax)
		{
			throw invalid_argument("Invalid range parameters. From get_grid_point3s(...)");
		}

		if (dlon <= 0 || dlat <= 0)
		{
			throw invalid_argument("Invalid interval parameters. From get_grid_point3s(...)");
		}

		int lonnum = floor((lonmax-lonmin)/dlon) + 1;
		int latnum = floor((latmax-latmin)/dlat) + 1;

		out_ps.resize(lonnum*latnum);
		for (int j = 0; j < latnum; j++)
		{
			for (int i = 0; i < lonnum; i++)
			{
				out_ps.at(j*lonnum+i).lon = lonmin + dlon*i;
				out_ps.at(j*lonnum+i).lat = latmin + dlat*j;
				out_ps.at(j*lonnum+i).rad = rad;
			}
		}
		return;
	}
}

#endif // _GCTL_POINT3S_H